#pragma once

#include "operators.hpp"

//
// alt ("algebraic template library")
//
//written by Ryo Hashimoto(Apr 17, 2008)

namespace alt{

    template<typename ValueType,size_t N>
    class algebraic_equation :
        public multiplicative<algebraic_equation<ValueType,N>,ValueType >,
        public additive<algebraic_equation<ValueType,N>,ValueType >
	{
    public:
        typedef ValueType value_type;
        typedef std::complex<value_type> solution_type;
        typedef std::vector<solution_type> solution_container;
        static const size_t order = N;
        typedef algebraic_equation<value_type,N> this_type;
    private:
        value_type m_a[N+1];
        
        //operator implementation detail
        template<bool this_type_is_result_type,size_t N_r>
        struct add_impl{
            typedef this_type lhs_type;
            typedef algebraic_equation<value_type,N_r> rhs_type;
            typedef rhs_type result_type;

            result_type operator()(const lhs_type& lhs,const rhs_type& rhs)const
            {
                result_type res(rhs);
                for(size_t i=0;i!=N+1;++i){
                    res[i] += lhs[i];
                }
                return res;
            }
        };
        template<size_t N_r>
        struct add_impl<true,N_r>{
            typedef this_type lhs_type;
            typedef algebraic_equation<value_type,N_r> rhs_type;
            typedef this_type result_type;

            result_type operator()(const lhs_type& lhs,const rhs_type& rhs)const
            {
                result_type res(lhs);
                for(size_t i=0;i!=N_r+1;++i){
                    res[i] += rhs[i];
                }
                return res;
            }
        };
        
        template<typename ResultType,typename EqType,size_t M>
        struct pow_impl{
            ResultType operator()(const EqType& eq)const
            {
                return eq*pow_impl<
                    algebraic_equation<typename EqType::value_type,EqType::order*(M-1)>,
                    EqType,
                    M-1
                >()(eq);
            }
        };
        template<typename ResultType,typename EqType>
        struct pow_impl<ResultType,EqType,0>{
            ResultType
                operator()(const EqType& eq)
            {
                ResultType res;
                res[0] = 1;
                return res;
            }
        };

    public:
        //constructor
        algebraic_equation()
        {
        }
        template<typename IterType>
        algebraic_equation(IterType begin,IterType end)
        {
            value_type* p = m_a;
            for(;begin!=end;++begin,++p){
                *p = *begin;
            }
        }

        //functional
        template<typename VariantType>
        solution_type operator()(const VariantType& x)const
        {
            solution_type res=0;

            for(size_t i=0;i!=N+1;++i){
                solution_type a = (*this)[i];
                for(size_t j=0;j!=i;++j){
                    a *= x;
                }
                res += a;
            }
            return res;
        }

        //accessor
        value_type& operator[](size_t index)
        {
            assert(index<N+1);
            return m_a[index];
        }
        const value_type& operator[](size_t index)const
        {
            assert(index<N+1);
            return m_a[index];
        }

        //add
        template<size_t M>
        algebraic_equation<value_type,((N>M)?N:M)>
            operator+(const algebraic_equation<value_type,M> &r)const
        {
            const size_t N_result = (N>M)?N:M;
            typedef algebraic_equation<value_type,N_result> result_type;
            return add_impl<N==N_result,M>()(*this,r);
        }
        this_type& operator+=(const value_type& r)
        {
            (*this)[0] += r;
            return *this;
        }

        //subtract
        template<size_t M>
        algebraic_equation<value_type,((N>M)?N:M)>
            operator-(const algebraic_equation<value_type,M> &r)const
        {
            return -r+(*this);
        }
        this_type& operator-=(const value_type& r)
        {
            (*this)[0] -= r;
            return *this;
        }
        this_type operator-()const
        {
            this_type res;

            for(size_t i=0;i!=N+1;++i){
                res[i] = -(*this)[i];
            }

            return res;
        }

        //multiply
        template<size_t M>
        algebraic_equation<value_type,N+M>
            operator*(const algebraic_equation<value_type,M> &r)const
        {
            const size_t N_result = N+M;
            typedef algebraic_equation<value_type,N_result> result_type;

            result_type res;

            for(size_t i=0;i!=N_result+1;++i)res[i]=0;

            for(size_t i=0;i!=N+1;++i){
                for(size_t j=0;j!=M+1;++j){
                    res[i+j] += (*this)[i]*r[j];
                }
            }

            return res;
        }
        this_type& operator*=(const value_type& r)
        {
            for(size_t i=0;i!=N+1;++i){
                (*this)[i] *= r;
            }
            return *this;
        }

        //divide
        this_type& operator/=(const value_type& r)
        {
            for(size_t i=0;i!=N+1;++i){
                (*this)[i] /= r;
            }
            return *this;
        }

        //equal
        template<size_t M>
        algebraic_equation<value_type,((N>M)?N:M)>
            operator==(const algebraic_equation<value_type,M> &r)const
        {
            const size_t N_result = (N>M)?N:M;
            typedef algebraic_equation<value_type,N_result> result_type;

            return (*this)-r;
        }
        this_type operator==(const value_type& r)
        {
            return *this-r;
        }

        //power
        template<size_t M>
        algebraic_equation<ValueType,N*M> pow()const
        {
            typedef algebraic_equation<ValueType,N*M> result_type;
            return pow_impl<result_type,this_type,M>()(*this);
        }
    };
}
